The PIES model
cp
Gurobi
AMPL
short
Ax = q(p),
s.t. Bx = b,
x >= 0
where the prices p are the duals to the first constraint.
MCP: A'p + B'v + c >= 0, x >= 0, comp.
-Ax + q(p) = 0, p free, comp.
-Bx + b = 0, v free, comp.
Of course, the variables x and v are going to be
split up further in the model.
References:
William W. Hogan, Energy Policy Models for Project Independence,
Computers \& Operations Research (2), 1975.
N. Josephy, A Newton Method for the PIES energy model,
Tech Report, Mathematics Research Center, UW-Madison, 1979.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */
set comod ordered; /* coal and light and heavy oil */
set R; /* resources */
data;
set comod := C L H;
set R := Capital Steel;
model;
set creg := 1 .. 2; /* coal producing regions */
set oreg := 1 .. 2; /* crude oil producing regions */
set ctyp := 1 .. 3; /* increments of coal production */
set otyp := 1 .. 2; /* increments of oil production */
set refin := 1 .. 2; /* refineries */
set users := 1 .. 2; /* consumption regions */
param rmax {R}; /* maximum resource usage */
param cmax {creg,ctyp}; /* coal prod. limits */
param omax {oreg,otyp}; /* oil prod. limits */
param rcost {refin}; /* refining cost */
param q0 {comod}; /* base demand for commodities */
param p0 {comod}; /* base prices for commodities */
param demand {comod,users}; /* computed at optimality */
param output {refin,comod}; /* % output of light/heavy oil */
param esub {comod,comod}; /* cross-elasticities of substitution */
param cruse {R,creg,ctyp}; /* resource use in coal prod */
param oruse {R,oreg,otyp}; /* resource use in oil prod */
param ccost {creg,ctyp}; /* coal prod. cost */
param ocost {oreg,otyp}; /* oil prod. cost */
param ctcost {creg,users};
param otcost {oreg,refin};
param ltcost {refin,users}; /* light oil trans costs */
param htcost {refin,users}; /* heavy oil trans costs */
param i_c {creg,ctyp};
param i_o {oreg,otyp};
param i_ct {creg,users}; /* initial trans */
param i_ot {oreg,refin}; /* initial trans */
param i_lt {refin,users}; /* initial trans */
param i_ht {refin,users}; /* initial trans */
param iprice {comod,users}; /* initial price estimate */
var c {cr in creg,ct in ctyp} := i_c[cr,ct]; /* coal production */
var o {r in oreg,ot in otyp} := i_o[r,ot]; /* oil production */
var ct {cr in creg,u in users} := i_ct[cr,u]; /* coal transportation levels */
var ot {Or in oreg,r in refin} := i_ot[Or,r]; /* crude oil transportation levels */
var lt {r in refin,u in users} := i_lt[r,u]; /* light transportation levels */
var ht {r in refin,u in users} := i_ht[r,u]; /* heavy transportation levels */
var p {co in comod,u in users} := iprice[co,u]; /* commodity prices */
var mu {R} := 1; /* dual to ruse cons.; marginal utility */
var cv {creg} := 1; /* dual to cmbal */
var ov {oreg} := 1;
var lv {refin} := 1;
var hv {refin} := 1;
delc {cr in creg,t in ctyp}:
0 <= c[cr,t] <= cmax[cr,t] complements
ccost[cr,t] + (sum {res in R} cruse[res,cr,t]*mu[res]) - cv[cr];
delo {r in oreg,t in otyp}:
0 <= o[r,t] <= omax[r,t] complements
ocost[r,t] + (sum {res in R} oruse[res,r,t]*mu[res]) - ov[r];
delct {cr in creg,u in users}:
0 <= ct[cr,u] complements
ctcost[cr,u] + cv[cr] >= p["C",u];
delot {Or in oreg,r in refin}:
0 <= ot[Or,r] complements
otcost[Or,r] + rcost[r] + ov[Or] >=
output[r,"L"] * lv[r] + output[r,"H"] * hv[r];
dellt {r in refin,u in users}:
0 <= lt[r,u] complements
ltcost[r,u] + lv[r] >= p["L",u];
delht {r in refin,u in users}:
0 <= ht[r,u] complements
htcost[r,u] + hv[r] >= p["H",u];
dembal {co in comod,u in users}: /* excess supply of product */
.1 <= p[co,u] complements
(if ord(co) = 1 then sum {r in creg} ct[r,u])
+ (if ord(co) = 2 then sum {r in refin} lt[r,u])
+ (if ord(co) = 3 then sum {r in refin} ht[r,u])
>= q0[co] * prod {cc in comod} (p[cc,u]/p0[cc])**esub[co,cc];
cmbal {cr in creg}: /* coal material balance */
cv[cr] complements
sum {t in ctyp} c[cr,t] = sum {u in users} ct[cr,u];
ombal {Or in oreg}: /* oil material balance */
ov[Or] complements
sum {t in otyp} o[Or,t] = sum {r in refin} ot[Or,r];
lmbal {r in refin}: /* light material balance */
lv[r] complements
sum {Or in oreg} ot[Or,r] * output[r,"L"] =
sum {u in users} lt[r,u];
hmbal {r in refin}: /* heavy material balance */
hv[r] complements
sum {Or in oreg} ot[Or,r] * output[r,"H"] =
sum {u in users} ht[r,u];
ruse {res in R}: /* resource use constraints */
0 <= mu[res] complements
rmax[res] >=
sum {cr in creg, t in ctyp} c[cr,t]*cruse[res,cr,t]
+ sum {Or in oreg, t in otyp} o[Or,t]*oruse[res,Or,t];
]]>
pies.dat
param rmax :=
Capital 35000
Steel 12000;
param cmax :=
1 1 300
1 2 300
1 3 400
2 1 200
2 2 300
2 3 600;
param omax:=
1 1 1100
1 2 1200
2 1 1300
2 2 1100;
param rcost :=
1 6.5
2 5;
param q0 :=
C 1000
L 1200
H 1000;
param p0 :=
C 12
L 16
H 12;
param output:=
1 L .6
1 H .4
2 L .5
2 H .5;
param esub
: C L H :=
C -.75 .1 .2
L .1 -.5 .2
H .2 .1 -.5 ;
param cruse
[Capital,*,*]
: 1 2 3 :=
1 1 5 10
2 1 5 6
[Steel,*,*]
: 1 2 3 :=
1 1 2 3
2 1 4 5 ;
param oruse
[Capital,*,*]
: 1 2 :=
1 0 10
2 0 15
[Steel,*,*]
: 1 2 :=
1 0 4
2 0 2 ;
param ccost
: 1 2 3 :=
1 5 6 8
2 4 5 7 ;
param ocost
: 1 2 :=
1 1 1.5
2 1.25 1.5 ;
param ctcost
: 1 2 :=
1 1 2.5
2 .75 2.75 ;
param otcost
: 1 2 :=
1 2 3
2 4 2 ;
param ltcost
: 1 2 :=
1 1 1.2
2 1 1.5 ;
param htcost
: 1 2 :=
1 1 1.2
2 1 1.5 ;
param i_c
: 1 2 3 :=
1 300 300 400
2 200 300 600 ;
param i_o :=
: 1 2 :=
1 1100 1000
2 1300 1000 ;
param i_ct
: 1 2 :=
1 0 828
2 1016 84 ;
param i_ot
: 1 2 :=
1 2075 0
2 0 2358 ;
param i_lt
: 1 2 :=
1 22 1223
2 1179 0 ;
param i_ht
: 1 2 :=
1 0 830
2 998 180 ;
param iprice
: 1 2 :=
C 11.7 13.7
L 15.8 16.0
H 11.9 12.4 ;
]]>
print " Solving the pies model:";
solve;
print " The commodity prices are:";
display p;